%%%%%%%%%%%%%%%%%%%SV_simulation%%%%%%%%%%%%
%Monte Carlo simulation for stochastic volatility model
%Original equations: 
% ds = mu*s dt + sigma1*sqrt(v)*s dw1;
% dv = k*(theta-v) dt + sigma2*sqrt(v) dw2; cov(w1,w2)=p;
%Euler dissociation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Computation of Gaussian variables
DeltaT=0.001;NbTraj=1;NbStep=1/DeltaT;
%NbTraj: Number of trajectories
%NbStep: Number of step to simulate
X0=1;a=0.03;k=10;b1=0.4;b2=0.8;
c=0.4;p=-0.7;v0=0.4;
%a=drift,b=sigma,c=theta,p=correlation coefficient
DeltaW1 = sqrt(DeltaT)*randn(NbTraj,NbStep);
DeltaW2 = sqrt(DeltaT)*randn(NbTraj,NbStep);
%v(t):Matrices used to store the approximations
vraj=zeros(NbTraj,NbStep+1);
vraj(:,1) = v0*ones(NbTraj,1);
%v(t):Loop that simulates all the time steps
for i=1:NbStep
    vraj(:,i+1) = vraj(:,i).*(1-k.*DeltaT)+b2*sqrt(vraj(:,i)).*(p.*DeltaW1(:,i)...
    +sqrt(1-p^2).*DeltaW2(:,i))+k*c.*DeltaT;
end
%s(t):Matrices used to store the approximations
Traj=zeros(NbTraj,NbStep+1);
Traj(:,1) = X0*ones(NbTraj,1);
%s(t):Loop that simulates all the time steps
for i=1:NbStep
    Traj(:,i+1) = Traj(:,i).*(1+a*DeltaT+b1*sqrt(vraj(:,i)).*DeltaW1(:,i));
end
%Draw the pic
subplot(1,2,1),plot(vraj);
xlabel('t');
ylabel('V(t)');
subplot(1,2,2),plot(Traj);
xlabel('t');
ylabel('S(t)');
set(gcf,'unit','centimeters','position',[10 5 20 13]);
